knitr::opts_chunk$set(echo = TRUE)
list.of.packages <- c("gsheet", "tidyverse", "janitor", "plotly", "glmmTMB", "metafor", "broom.mixed", "car", "ggpubr", "scales") #add new libraries here
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# Load all libraries
lapply(list.of.packages, FUN = function(X) {
do.call("require", list(X))
})
sessionInfo()
NOTE: data is read directly from the GoogleSheet using a share link that was set to “anyone with a link can view”
# Read in data from GoogleSheet
data.fert <- as_tibble(gsheet2tbl('https://docs.google.com/spreadsheets/d/111SuH548Et6HDckjbQjtYk-B8A5VfNkUYZgcUNRPxxY/edit?usp=sharing'))
## Warning: Missing column names filled in: 'X22' [22]
#replace NR (not reported) with NA and convert columns to factor / numeric where needed
data.fert <- data.fert %>%
na_if("NR") %>%
clean_names() %>% # fill in spaces with underscores for column names
mutate_at(c('phylum', 'study', 'taxonomic_group', 'common_name', 'latin_name', 'error_statistic'), as.factor) %>%
mutate_at(c('p_h_experimental', 'p_h_control', 'p_co2_experimental', 'p_co2_control', 'ave_fert_percent_p_h', 'error_percent_p_h', 'number_trials_p_h', 'insemination_time', 'sperm_per_m_l', 'sperm_egg_ratio', 'number_females', 'number_males'), as.numeric)
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
#data.fert %>% select(study, ave_fert_percent_p_h, error_percent_p_h, error_statistic, number_trials_p_h) %>% View()
ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) +
geom_point(size=1.5, width=0.02) +
facet_wrap(~phylum, scale="free") +
geom_smooth(method="lm", se=TRUE, aes(fill=taxonomic_group)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum, linear") +
theme_minimal()
)
## Warning: Ignoring unknown parameters: width
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=study, col=study, text=`common_name`)) +
geom_point(size=1.5, width=0.02) +
facet_wrap(~phylum, scale="free") +
geom_smooth(method="lm", se=TRUE, aes(fill=taxonomic_group)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum, linear") +
theme_minimal()
)
## Warning: Ignoring unknown parameters: width
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
## Warning in qt((1 - level)/2, df): NaNs produced
ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) +
geom_point(size=1, width=0.02) +
facet_wrap(~phylum, scale="free") +
geom_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 2, raw=TRUE), aes(fill=taxonomic_group)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum, polynomial") +
theme_minimal()
)
## Warning: Ignoring unknown parameters: width
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
Upper 95% CI = Mean + 1.96*SE; I recorded the difference between the upper 95%CI and the mean (Upper 95%CI - Mean). To convert I will use:
SE = (Upper 95%CI - Mean) / 1.96
SD = ((Upper 95%CI - Mean) / 1.96) * sqrt(n)
SE= SD/sqrt(n)
SD = SE*sqrt(n)
where n=sample size
data.fert <- data.fert %>%
mutate(SE = case_when(error_statistic == "SD" ~ error_percent_p_h/sqrt(number_trials_p_h),
error_statistic == "95% CI" ~ error_percent_p_h/1.96,
is.na(error_statistic) ~ error_percent_p_h,
error_statistic == "SE" ~ error_percent_p_h))
data.fert <- data.fert %>%
mutate(SD = case_when(error_statistic == "SE" ~ error_percent_p_h*sqrt(number_trials_p_h),
error_statistic == "95% CI" ~ (error_percent_p_h/1.96)*sqrt(number_trials_p_h),
is.na(error_statistic) ~ error_percent_p_h,
error_statistic == "SD" ~ error_percent_p_h))
data.fert <- data.fert %>%
mutate(pH_delta = p_h_control-p_h_experimental)
data.fert <- data.fert %>%
mutate(ave_fert_proport = case_when(ave_fert_percent_p_h <= 100 ~ ave_fert_percent_p_h/100,
ave_fert_percent_p_h > 100 ~ 1))
data.fert$ave_fert_proport %>% hist()
ggplot(data.fert, aes(group=phylum, col=phylum)) + geom_density(aes(ave_fert_proport))
## Warning: Removed 13 rows containing non-finite values (stat_density).
# How many studies per ~phylum?
data.fert %>%
select(phylum, study) %>%
distinct(phylum, study) %>%
group_by(phylum) %>% count()
## # A tibble: 4 x 2
## # Groups: phylum [4]
## phylum n
## <fct> <int>
## 1 Cnidaria 6
## 2 Crustacea 5
## 3 Echinodermata 12
## 4 Mollusca 18
# How many studies per taxonomic group?
data.fert %>%
select(phylum, study, taxonomic_group) %>%
distinct(phylum, study, taxonomic_group) %>%
group_by(taxonomic_group) %>% count()
## # A tibble: 10 x 2
## # Groups: taxonomic_group [10]
## taxonomic_group n
## <fct> <int>
## 1 abalone 2
## 2 clam 5
## 3 copepod 3
## 4 coral 6
## 5 mussel 4
## 6 non-copepod 2
## 7 non-urchin 2
## 8 oyster 5
## 9 scallop 3
## 10 sea urchin 11
weights <- metafor::escalc(measure='MN',
mi=data.fert$ave_fert_percent_p_h,
sdi = data.fert$SD,
ni=data.fert$number_trials_p_h, options(na.action="na.pass"))
data.fert$w <-weights$vi
# How many studies per ~phylum?
data.fert %>% drop_na(w) %>%
select(phylum, study) %>%
distinct(phylum, study) %>%
group_by(phylum) %>% count()
## # A tibble: 4 x 2
## # Groups: phylum [4]
## phylum n
## <fct> <int>
## 1 Cnidaria 5
## 2 Crustacea 5
## 3 Echinodermata 11
## 4 Mollusca 15
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #yes
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 9.1418 1 0.002498 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 0.0811 1 0.7758
Anova(glmmTMB(ave_fert_proport ~ phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## phylum 1.9538 3 0.582
Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 0.0557 1 0.8134
Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_per_m_l 0.024 1 0.877
Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_egg_ratio 2e-04 1 0.99
Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_females 0.054 1 0.8162
Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_males 0.2089 1 0.6476
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*p_h_control + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 9.1034 1 0.002551 **
## p_h_control 0.0432 1 0.835310
## p_h_experimental:p_h_control 0.0804 1 0.776795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*insemination_time + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 10.6640 1 0.001092 **
## insemination_time 0.7586 1 0.383757
## p_h_experimental:insemination_time 0.8087 1 0.368510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 7.9631 1 0.004774 **
## phylum 4.2472 3 0.235975
## p_h_experimental:phylum 2.9060 3 0.406354
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*sperm_per_m_l + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 1.1881 1 0.2757
## sperm_per_m_l 0.0272 1 0.8691
## p_h_experimental:sperm_per_m_l 0.0009 1 0.9762
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*sperm_egg_ratio + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 3.6744 1 0.05526 .
## sperm_egg_ratio 0.0047 1 0.94507
## p_h_experimental:sperm_egg_ratio 0.0000 1 0.99590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*number_females + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 8.7577 1 0.003083 **
## number_females 0.0718 1 0.788736
## p_h_experimental:number_females 0.0786 1 0.779138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental*number_males + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 8.8671 1 0.002903 **
## number_males 0.1893 1 0.663490
## p_h_experimental:number_males 0.0179 1 0.893581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
car::Anova(best <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #phylum, pH, & phylum:pH sign. factors
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 9.1418 1 0.002498 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(best)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: drop_na(data.fert, w)
## Weights: 1/(w + 1)^2
##
## AIC BIC logLik deviance df.resid
## 73.7 84.8 -33.9 67.7 295
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 1.929 1.389
## Number of obs: 298, groups: study, 33
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.382 7.673 -2.917 0.00354 **
## p_h_experimental 3.098 1.025 3.023 0.00250 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.fert <- function(pH, model) {
linear.predictor <- as.vector((summary(best)$coefficients$cond[,"Estimate"]["(Intercept)"] +
summary(best)$coefficients$cond[,"Estimate"]["p_h_experimental"]*pH))
predicted <- exp(linear.predictor) / (1+exp(linear.predictor))
return(paste("Fertilization rate predicted for pH ", pH, ": ",
scales::percent(x=predicted, accuracy = .01), sep=""))
}
# Estimate % fertilization @ pH 8.0 using hand-typed equation
exp(-22.382104 + 3.098070*8) / (1 + exp(-22.382104 + 3.098070*8))
## [1] 0.9170144
# now use function
predict.fert(8.0, best)
## [1] "Fertilization rate predicted for pH 8: 91.70%"
predict.fert(7.5, best)
## [1] "Fertilization rate predicted for pH 7.5: 70.13%"
predict.fert(7.0, best)
## [1] "Fertilization rate predicted for pH 7: 33.28%"
predict.fert(6.0, best)
## [1] "Fertilization rate predicted for pH 6: 2.20%"
confint(best)
## 2.5 % 97.5 % Estimate
## cond.(Intercept) -37.4217617 -7.342447 -22.382104
## cond.p_h_experimental 1.0897904 5.106349 3.098070
## study.cond.Std.Dev.(Intercept) 0.5899624 3.270399 1.389033
aa5 <- augment(best, data=drop_na(data.fert, w))
## Warning in mr - fitted(object): longer object length is not a multiple of
## shorter object length
#ggplotly(
ggplot(aa5, aes(x=.fitted,y=.resid)) +
geom_point() +
geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
#)
ph.min.max <- drop_na(data.fert, w) %>%
select(phylum, p_h_experimental) %>%
group_by(phylum) %>%
summarize(min=min(p_h_experimental, na.rm=TRUE), max=max(p_h_experimental, na.rm=TRUE))
phylum.list <- list()
for (i in 1:nrow(ph.min.max)) {
phylum.list[[i]] <- data.frame(ph=c(seq(from=as.numeric(ph.min.max[i,"min"]),
to=as.numeric(ph.min.max[i,"max"]),
by=0.01)),
phylum=rep(c(ph.min.max[i,"phylum"])))
}
new.data <- bind_rows(phylum.list) %>% purrr::set_names(c("p_h_experimental", "phylum"))
new.data$study <- NA
new.data$w <- NA
predict.test.df <- predict(best, newdata = new.data, se.fit = TRUE, type="response")
predict.test.df.df <- predict.test.df %>%
as.data.frame() %>%
cbind(new.data)
#scales::show_col(c("#e41a1c","#4daf4a","#ff7f00","#984ea3",'#377eb8'))
Fertilization success (%) by experimental pH across marine taxa examined in this review. Meta-analysis was performed using a binomial regression model, and indicates that fertilization success decreases with pH across Crustacea (5 studies), Echinodermata (12 studies), and Mollusca (18 studies). Fertilization success was not significantly affected by pH in Cnidaria (4 studies). Each point reflects the average % fertilization reported by one study at an experimental pH.
# Examine pH-experimental significance as predictor using X-squared test
Anova(best)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 9.1418 1 0.002498 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#ggplotly(
ggplot() +
geom_jitter(data=drop_na(data.fert, w), aes(x=p_h_experimental, y=ave_fert_proport, group=phylum, col=phylum, label=study), size=1.2, width=0.03) + #, col="gray40"
#facet_wrap(~phylum, scales="free") +
theme_minimal() +
ggtitle("Fertilization success ~ pH with binomial-regression model predictions (all phyla)") +
xlab("Experimental pH") + ylab("Proportion fertilization success") +
scale_color_manual(name=NULL, values=c("#ca0020","#f4a582","#92c5de",'#0571b0')) +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.test.df.df, aes(x=p_h_experimental, y=fit), col="gray50") + #, col=phylum
geom_ribbon(data = predict.test.df.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="gray50") + theme(legend.position = "top") + guides(colour = guide_legend(override.aes = list(size=4))) +
annotate(geom="text", x=6.4, y=0.9, size=5, colour="gray50", label=paste("χ2 p-value =", round(Anova(best)[,"Pr(>Chisq)"], digits =4), sep=" ")) #) #add fill=phylum in geom_ribbon aes if color desired
## Warning: Ignoring unknown aesthetics: label
## Warning: Removed 12 rows containing missing values (geom_point).
show_col is from scales library)show_col(c("#ca0020","#f4a582","#92c5de",'#0571b0'))
#e41a1c = Cnidaria (red)
#ff7f00 = Crustacea (orange)
#4daf4a = Echinodermata (green)
#377eb8 = Mollusca (blue)
#ca0020 = coral
#f4a582 = crustacean
#92c5de = echinoderm
#0571b0 = mollusc
ggplotly(
ggplot() +
geom_jitter(data=drop_na(data.fert, w), aes(x=p_h_experimental, y=ave_fert_proport, group=phylum, col=phylum, label=study), size=1.2, width=0.03, col="gray40") +
facet_wrap(~phylum, scales="free") + theme_minimal() +
ggtitle("Fertilization success ~ pH with binomial-regression model predictions") +
xlab("Experimental pH") + ylab("Proportion fertilization success") +
#scale_color_manual(name=NULL, values=c("#e41a1c","#ff7f00","#4daf4a",'#377eb8')) +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.test.df.df, aes(x=p_h_experimental, y=fit), col="gray50") + #, col=phylum
geom_ribbon(data = predict.test.df.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="gray50") + theme(legend.position = "none")) #add fill=phylum in geom_ribbon aes if color desired
## Warning: Ignoring unknown aesthetics: label
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 1.5152 1 0.2183
Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 1.555 1 0.2124
Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 0.1092 1 0.7411
Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_per_m_l 1
Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_egg_ratio 1
Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_females 4e-04 1 0.9844
Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_males 0.0058 1 0.9392
# None significant, but still develop model for plot
model.cnidaria <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.cnidaria)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 1.5152 1 0.2183
summary(model.cnidaria)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Cnidaria")
## Weights: 1/(w + 1)^2
##
## AIC BIC logLik deviance df.resid
## 28.1 34.0 -11.1 22.1 50
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.4241 0.6513
## Number of obs: 53, groups: study, 5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -41.697 34.526 -1.208 0.227
## p_h_experimental 5.406 4.392 1.231 0.218
paste("Experimental pH χ2 p-value =", round(Anova(model.cnidaria)[,"Pr(>Chisq)"], digits = 3), sep=" ")
## [1] "Experimental pH χ2 p-value = 0.218"
predict.cnidaria <- predict(model.cnidaria, newdata = subset(new.data, phylum=="Cnidaria"), se.fit = TRUE, type="response")
predict.cnidaria.df <- predict.cnidaria %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Cnidaria"))
plot.cnidaria <- ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>%
filter(phylum=="Cnidaria"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#ca0020") +
ggtitle("Cnidaria") +
xlab(NULL) + ylab("Proportion fertilization success") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.cnidaria.df, aes(x=p_h_experimental, y=fit), col="#ca0020") +
geom_ribbon(data = predict.cnidaria.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#ca0020") +
annotate(geom="text", x=6.5, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.cnidaria)[,"Pr(>Chisq)"], digits = 3), sep=" ")) #Significance of experimental pH \nas predictor:
print(plot.cnidaria)
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 0.0635 1 0.801
Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; function evaluation
## limit reached without convergence (9). See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 0.0028 1 0.9579
Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; extreme or very small
## eigen values detected. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 1
model.crustacea <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
Anova(model.crustacea)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 0.0635 1 0.801
summary(model.crustacea)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Crustacea")
## Weights: 1/(w + 1)^2
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 21
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 7.022e-109 8.38e-55
## Number of obs: 24, groups: study, 5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.51 75.22 -0.273 0.785
## p_h_experimental 2.58 10.24 0.252 0.801
predict.crustacea <- predict(model.crustacea, newdata = subset(new.data, phylum=="Crustacea"), se.fit = TRUE, type="response")
predict.crustacea.df <- predict.crustacea %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Crustacea"))
plot.crustacea <- ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>%
filter(phylum=="Crustacea"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#f4a582") +
ggtitle("Crustacea") +
xlab(NULL) + ylab(NULL) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.crustacea.df, aes(x=p_h_experimental, y=fit), col="#f4a582") +
geom_ribbon(data = predict.crustacea.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#f4a582") +
annotate(geom="text", x=6.5, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.crustacea)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.crustacea)
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #yes
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 6.3225 1 0.01192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 0.2848 1 0.5936
Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 3.7605 1 0.05248 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_per_m_l 0 1 0.9963
Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_egg_ratio 0.0671 1 0.7957
Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_females 0.2046 1 0.651
Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_males 2.4125 1 0.1204
model.echinodermata <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.echinodermata)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 6.3225 1 0.01192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.echinodermata)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Echinodermata")
## Weights: 1/(w + 1)^2
##
## AIC BIC logLik deviance df.resid
## 29.0 37.9 -11.5 23.0 142
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 2.661 1.631
## Number of obs: 145, groups: study, 11
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -41.938 16.894 -2.482 0.0131 *
## p_h_experimental 5.619 2.235 2.514 0.0119 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.echinodermata <- predict(model.echinodermata, newdata = subset(new.data, phylum=="Echinodermata"), se.fit = TRUE, type="response")
predict.echinodermata.df <- predict.echinodermata %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Echinodermata"))
plot.echino <- ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>%
filter(phylum=="Echinodermata"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#92c5de") +
ggtitle("Echinodermata") +
xlab("Experimental pH") + ylab("Proportion fertilization success") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.echinodermata.df, aes(x=p_h_experimental, y=fit), col="#92c5de") +
geom_ribbon(data = predict.echinodermata.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#92c5de") +
annotate(geom="text", x=6.4, y=0.85, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.echinodermata)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.echino)
Anova(glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #meh
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 2.8595 1 0.09084 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(glmmTMB(ave_fert_proport ~ p_h_control + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
## NaN function evaluation
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 1
Anova(glmmTMB(ave_fert_proport ~ insemination_time + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 0.008 1 0.9288
Anova(glmmTMB(ave_fert_proport ~ sperm_per_m_l + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_per_m_l 1
Anova(glmmTMB(ave_fert_proport ~ sperm_egg_ratio + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in eval(family$initialize): Model convergence problem; non-positive-
## definite Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_egg_ratio 1
Anova(glmmTMB(ave_fert_proport ~ number_females + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_females 0.5883 1 0.4431
Anova(glmmTMB(ave_fert_proport ~ number_males + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_males 0.3482 1 0.5552
model.mollusca <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.mollusca)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 2.8595 1 0.09084 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.mollusca)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Mollusca")
## Weights: 1/(w + 1)^2
##
## AIC BIC logLik deviance df.resid
## 19.8 26.8 -6.9 13.8 73
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 2.209e-06 0.001486
## Number of obs: 76, groups: study, 14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.5492 6.4052 -1.491 0.1360
## p_h_experimental 1.5196 0.8986 1.691 0.0908 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.mollusca <- predict(model.mollusca, newdata = subset(new.data, phylum=="Mollusca"), se.fit = TRUE, type="response")
predict.mollusca.df <- predict.mollusca %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Mollusca"))
plot.mollusca <- ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>%
filter(phylum=="Mollusca"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#0571b0") +
ggtitle("Mollusca") +
xlab("Experimental pH") + ylab(NULL) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, size=13, colour="gray30")) +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.mollusca.df, aes(x=p_h_experimental, y=fit), col="#0571b0") +
geom_ribbon(data = predict.mollusca.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#0571b0") +
annotate(geom="text", x=6.4, y=0.9, size=3.5, colour="gray20", label=paste("χ2 p-value =\n", round(Anova(model.mollusca)[,"Pr(>Chisq)"], digits = 3), sep=" "))
print(plot.mollusca)
## Warning: Removed 12 rows containing missing values (geom_point).
all.plots <- ggarrange(plot.cnidaria + ylab(NULL), plot.crustacea, plot.echino + ylab(NULL), plot.mollusca + rremove("y.text"), ncol=2, nrow=2)
## Warning: Removed 12 rows containing missing values (geom_point).
annotate_figure(all.plots, left = text_grob("Proportion Fertilization Success", color = "gray20", rot = 90))
#add the following to remove x and y axis labels: ` + rremove("x.text") + rremove("y.text")`
# add the following to add plot labels: `, labels=c("A", "B", "C", "D")`
data.fert <- data.fert %>%
mutate_at(vars(number_trials_p_h, insemination_time, sperm_per_m_l, sperm_egg_ratio, number_females, number_males), as.numeric)
car::Anova(covariates1 <- glmmTMB(ave_fert_proport ~ p_h_control*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no interaction
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_control 0.0001 1 0.992776
## p_h_experimental 9.4009 1 0.002169 **
## p_h_control:p_h_experimental 0.7993 1 0.371315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplotly(
ggplot(data.fert, aes(x=p_h_control, y=ave_fert_proport)) +
geom_jitter(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + ggtitle("% Fertilization ~ Control pH\nsign. interaction") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))
car::Anova(covariates2 <- glmmTMB(ave_fert_proport ~ insemination_time*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #no interaction
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## insemination_time 0.7835 1 0.3760640
## p_h_experimental 10.8791 1 0.0009725 ***
## insemination_time:p_h_experimental 0.7018 1 0.4021760
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplotly(
ggplot(data.fert, aes(x=insemination_time, y=ave_fert_proport)) +
geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) + ggtitle("% Fertilization ~ Insemination Time\nsign. interaction_") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))
car::Anova(covariates3 <- glmmTMB(ave_fert_proport ~ sperm_per_m_l*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_per_m_l 0.0121 1 0.9124
## p_h_experimental 1.3803 1 0.2400
## sperm_per_m_l:p_h_experimental 0.0004 1 0.9849
ggplotly(ggplot(data.fert, aes(x=sperm_per_m_l, y=ave_fert_proport)) +
geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) +
ggtitle("% Fertilization ~ Sperm Concentration\nnot sign.") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))
car::Anova(covariates4 <- glmmTMB(ave_fert_proport ~ sperm_egg_ratio*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## sperm_egg_ratio 0.0073 1 0.931826
## p_h_experimental 6.8239 1 0.008994 **
## sperm_egg_ratio:p_h_experimental 0.0000 1 0.995566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplotly(ggplot(data.fert, aes(x=sperm_egg_ratio, y=ave_fert_proport)) +
geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) +
ggtitle("% Fertilization ~ Sperm:Egg Ratio\nnot sign.") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))
car::Anova(covariates5 <- glmmTMB(ave_fert_proport ~ number_females*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign.
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_females 0.0001 1 0.991430
## p_h_experimental 9.2586 1 0.002344 **
## number_females:p_h_experimental 0.0550 1 0.814625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplotly(ggplot(data.fert, aes(x=number_females, y=ave_fert_proport)) +
geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) +
ggtitle("% Fertilization ~ No. Females\nno interaction") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))
car::Anova(covariates6 <- glmmTMB(ave_fert_proport ~ number_males*p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = 1/(w+1)^2)) #not sign.
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## number_males 0.0535 1 0.817025
## p_h_experimental 9.2455 1 0.002361 **
## number_males:p_h_experimental 0.0423 1 0.836968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplotly(ggplot(data.fert, aes(x=number_males, y=ave_fert_proport)) +
geom_point(aes(color = p_h_experimental, label=common_name), size=2.2, pch=19) +
ggtitle("% Fertilization ~ No. Males\nsign. interaction") +
scale_color_gradientn(colours = rainbow(5)))
## Warning: Ignoring unknown aesthetics: label
# scale_color_gradient(low = "red", high = "blue"))